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We explore the role of non-Gaussian fluctuations in primordial black hole (PBH) formation 
and show that the standard Gaussian assumption, used in all PBH formation papers to date, is not 
justified. Since large spikes in power are usually associated with flat regions of the inflaton potential, 
quantum fluctuations become more important in the field dynamics, leading to mode-mode coupling 
and non-Gaussian statistics. Moreover, PBH production requires several sigma (rare) fluctuations 
in order to prevent premature matter dominance of the universe, so we are necessarily concerned 
with distribution tails, where any intrinsic skewness will be especially important. We quantify this 
argument by using the stochastic slow-roll equation and a relatively simple analytic method to obtain 
the final distribution of fluctuations. We work out several examples with toy models that produce 
PBH's, and test the results with numerical simulations. Our examples show that the naive Gaussian 
assumption can result in errors of many orders of magnitude. For models with spikes in power, 
our calculations give sharp cut-offs in the probability of large positive fluctuations, meaning that 
Gaussian distributions would vastly over-produce PBH's. The standard results that link inflation- 
produced power spectra and PBH number densities must then be reconsidered, since they rely 
quite heavily on the Gaussian assumption. We point out that since the probability distributions 
depend strongly on the nature of the potential, it is impossible to obtain results for general models. 
However, calculating the distribution of fluctuations for any specific model seems to be relatively 
straightforward, at least in the single inflaton case. 



I. INTRODUCTION 



Primordial black holes (PBH's) represent a link between quantum field theory and general relativity, both in 
their birth in inflation and their decay via Hawking radiation [Q. Beyond the exciting physics of their formation 
and evolution, PBH's have the status of a dark matter (DM) candidate. There is, in fact, renewed interest in this 
possibility, as Jedamzik j|] has recently investigated whether PBH formation during the QCD epoch could explain 
the existence of the ~ 0.5M Q MACHO's indicated by recent observations Q. Through the requirement that they not 
be over-produced, PBH's are also important as one of the few constraints on the inflationary Universe scenario 
Since these objects are an important potential link to the early universe, their formation should be examined in 
detail, especially if the standard Gaussian assumption about their formation may be producing errors of many orders 
in magnitude. 

In order to determine the degree to which non-Gaussian fluctuations are important, we use a stochastic inflation 
calculation to investigate several PBH-producing toy models. We solve for the probability distributions both ana- 
lytically and numerically, and show that the Gaussian assumption is significantly flawed. Specifically, for models 
associated with spikes in small-scale power, the Gaussian assumption over-estimates PBH production by many orders 
of magnitude. These results suggest that if some model of inflation is in danger of being ruled out due to overpro- 
duction of PBH's (e.g., as considered by (111), then looking closely for non-Gaussian behavior may serve to save the 
model, or at least alter the ranges of viable parameters. 

We start with a heuristic discussion. Our intuition about non-Gaussian statistics for single inflaton models is fostered 
by the simple relation between the inflation-produced power spectrum and the field dynamics driving inflation. Using 
this relation, we are able to make general statements about the behavior of inflaton fluctuations in models which form 
PBH's, and see qualitatively why non-Gaussian statistics might be important. In order to make this point, we briefly 
review the standard inflation scenario, set our notation, and use these results to illustrate why non-Gaussianity seems 
likely in PBH formation. 
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The inflationary Universe scenario provides the seeds of structure formation by linking initial density perturba- 
tions to quantum fluctuations in one or more scalar fields. See B for a review. During inflation, the energy density of 
the Universe, p, becomes dominated by a scalar field potential, V(4>), and the scale factor R expands superluminally 
(R ~ t n , n > 1). So length scales of fluctuations grow more quickly than the horizon, eventually passing out of causal 
contact, and only cross back inside the horizon after the inflation epoch has ended. In the chaotic scenario M, the 
spatially homogeneous field, <p, is initially displaced from the minimum of its potential and rolls downward with its 
motion governed by the Klein-Gordon and Einstein equations 

^ + 3H(4>)4> = -v'{4>), (l) 



H2 =3^r^ (vw+m ' (2) 

where m p i is the Planck mass and H is the Hubble parameter. In most instances, H sw V{4>) 1 ^ 2 and the 

"slow-roll" conditions 
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apply. The evolution of (f> is then friction dominated: 



V 
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and the universe continues to inflate until either of the slow roll conditions breaks down at some value <f> = 4> en d- 

As the physical size of fluctuations grows and they cross outside the scale of the horizon during inflation, their 
amplitude is set by quantum fluctuations in tf>. Since the energy density p is dominated by V(<f>), we have Sp ~ 
V'S(f> ~ V'H. After inflation ends R ~ t n ,n < 1 and perturbations begin to cross back inside the horizon. The rms 
magnitude of dp/ p when a scale re-enters the horizon, 5h, is simply related to its properties as it left the horizon. 
Using the gauge- invariant variable £ ~ Sp/(p + p), we have Q 
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where L is the co-moving scale of horizon crossing and the right hand side is evaluated when the fluctuation left the 
horizon at <fi = <fi ou t(L). When a fluctuation crosses outside the horizon it has physical size H~ 1 (4> out ), so to find the 
co-moving size today, we must scale by the amount of expansion since that time: L = [R /R(<fi 0U t)]H (<f> ou t)- We 
find the relation between 4> ut & n d the present-day length scale L by first determining the amount of expansion that 
occurred during the rest of inflation [R{<pend) / R{4>out)\ and then using entropy conservation to determine the amount 
of expansion that took place up to now [R / R(<fiend)]- We know R/R = H(<p) d In R — Hd(j>/<j) during inflation, so 
the inflationary expansion was quasi-exponential: 

£Ww) = f^M (-jSjr^) H-\<p out ) = ^exp(iV(^ Mt ))i?- 1 (^) (5) 

\ -fiend / out) J -fiend 

where the number of "e- folds" between exiting the horizon and the end of inflation, N(<fi ou t), is given by 

H{<j))d<j) 8tt Z^"" V(0)d(j) 



N(4> ovt ) = / « — / -^-f . (6) 

For the expansion since the end of inflation, we will here assume a "quick re-heat" approximation and use R /R en d ~ 
[3.17](T end /T ) w U(</) end ) 1 / 4 /0.85-ft:, where the numerical term arises due to the effective number of relativistic 
species. 



1 The expression for 8h is correct up to factors of order unity which depend e.g. on whether there is radiation or matter 
domination when the scale re-enters the horizon. 
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The crucial point from the above discussion is that every region of the inflaton potential is mapped directly to some 
scale of the fluctuation spectrum using equations (4)-(6). For typical model parameters, L ~ 1 Mpc corresponds to 
N((f>) ss 51 and L ~ 1 pc corresponds to AT(</>) w 37. We invert such relations to determine <^>(pc), <^>(Mpc), etc. 

Specific characteristics of the density spectrum over a range of wavelengths then allow us to make general statements 
about the behavior of the potential in some region of cf> space. As we shall see in Sec. II, PBH production requires 
a specially-tuned fluctuation spectrum: the fluctuation amplitude on small scales must increase over its value on 
the COBE scale by a factor of ~ 10 3 in order to produce an appreciable fraction of PBH's. It is this that makes 
non-Gaussianity important in these models. 

Hodges et al. |10| have examined the question of non-Gaussianity in some detail^] and showed that non-Gaussianity 



is negligible in single-inflaton models with spectra meeting the large-scale structure requirement Sh ~ 3 x 10~ 5 . 
However, models that produce PBH's must have density Sh> 0.01 over some range of small-scale wavelengths P,p"2|. 
In order to get a qualitative feel for why non-Gaussian fluctuations may be important in these models, consider a 
background perturbation expanded in terms of the usual creation and annihilation operators a\,, a k 

4>(x, t) = cj> c i{t) + 8cj>{x, t) 

5<p(x,t) = - r -^r { d 3 k [a k S<f> k e lkx + h.c] . 



(27T) 

Inserting the perturbed field into the free-field Klein-Gordon equation we obtain 

k 2 

S<Pk ~ jpHk + 3H5cf> k + V'Sfa = 0, (7) 

where we have expanded to first order in S<p and matched small terms. Gaussian statistics are preserved to first order 
since in the linear approximation all Fourier modes remain separate. 

As 5<ft becomes larger, however, the linear approximation breaks down and terms of order S(j) 2 become important. 
Any non-linearity implies mode-mode coupling and introduces non-Gaussianity. So the possibility of a non-Gaussian 
distribution of fluctuations increases with Scj). Moreover, PBH production requires fluctuations of several a above 
8 rms = Sh in order to keep their production relatively rare (see Sec. II for a discussion). So we may expect non- 
Gaussian statistics to play a role in PBH production in inflation. 

A much stronger effect becomes apparent if we examine the general nature of a potential V((f) that produces 
PBH's. Since PBH production requires a large increase in power on small scales, expression (4) requires that </> (or 
similarly V') become extremely small over some region of the potential. In these cases, writing down a decomposition 
like equation (7) makes no sense: every term is "small," so matching terms of equivalent size does not simplify the 
equation. In other words, quantum fluctuations become important in the dynamics of cj). We must keep 5<f> linked to 
<p, and in the slow-roll approximation we have 

The evolution becomes extremely non-linear and non-Gaussian effects become important. 

In what follows, we investigate the question of non-Gaussian statistics and PBH production in detail. In Sec. II 
we present the basics of PBH creation and the implications for the inflation-produced power spectrum. We review 
stochastic inflation in Sec. Ill and present the method of deriving probability distributions that we will use in our 
examples. In Sec. IV, we present three examples of inflationary scenarios which create PBH's in significant numbers, 
and solve for the distribution of fluctuations both analytically and numerically. We reserve Sec. V for conclusions 
and speculation. 

II. PRIMORDIAL BLACK HOLES AND CONSTRAINTS ON INFLATION 

A. Formation 

A PBH forms when a collapsing over-dense region is large enough to overcome the opposing force of pressure and 
falls within its Schwarzschild radius p^ ]. We review the basics of this process for a radiation-dominated universe, 



2 For examinations of specific models using the stochastic inflation approach see e.g. [ p"l| . 



3 



where the sound speed is given by c 8 = c/y/3, and the pressure is p — c 2 p = p/3. Consider a spherically symmetric^] 
region with density p greater than the background density p. The high-density region is governed by the positive 
curvature (K > 0) Friedmann equation 



H\t) = ^pii) - * (8) 



while the background space evolves as 



pi 



Following the standard approach, choose initial coordinates t = t = ij such that the initial expansion rates in each 
region are equal H = H = Hi. Using (8) and (9) we see that the initial size of the perturbation Si in terms of 
R(ti) = Ri is given by 



1 



Pl H 2 R 2 ' 



(10) 



Before the over-dense region breaks away from the expansion, R-i 1 / 2 ||, and we can estimate the time, t c , when 
the perturbed region stops expanding (H(t c ) = 0): 

t c ^ s-H t (11) 

and the corresponding expansion factor at this time 

R(t c ) = R C ~ S; 1/2 R,. (12) 

The perturbed region continues to contract only if it contains enough matter to overcome the pressure, R c > 
Rjeans = c s t c ~ "^^c- With this condition met, a black hole will inevitably form: we need not worry about vorticity 
or turbulence interfering with the formation process since the region will fall inside its Schwarzschild radius soon after 
turn-around (Rs = 2GM ~ 2Gp c R\ ~ Rlt~ 2 > \R C ). The requirement for PBH formation is then 

-U c < R c < t c (13) 
v3 

where the upper bound prevents the region of collapse from being larger than the horizon scale, leading to the 
formation of a separate universe [l^| . Now, dividing (13) by t c , we obtain a condition on R c /t c ~ Ri/ti ~ S^ 2 . This 
expression scales like a constant with time, so it is convenient to evaluate it at horizon crossing, R = t. We then have 
our final condition for perturbations at horizon crossing 

^< <*< 1- (14) 

Thus, given an initial fluctuation with 8 satisfying (14) at horizon crossing, we expect to form a PBH with mass M 
corresponding roughly to the horizon mass at that time (numerical results agree Jl5[| ) . 

B. Probing the Power Spectrum 

Although quite elegant in solving the problems of the standard big bang, inflation provides us only with a framework 
of ideas and not with any exact predictions for universal evolution. Depending on the scalar potential, number of 
fields, etc., inflation can produce a wide variety of fluctuation spectra e.g. ]To| , ^6| [l8[ and even the possibility of a 
low il Universe p9|. Because of the freedom in inflationary models, any means of constraining a particular scenario 



3 Note that the spherical assumption is well justified |l3) for many-cr "rare" fluctuations in a Gaussian random field. This 
result may actually be incorrect for non-Gaussian fluctuations and should be explored further. 
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becomes extremely important. The Cosmic Background Explorer (COBE) measurements and large-scale structure 
observations indicate Sh ~ 3 x 1CP 5 on large scales (~ 10 3 Mpc). These results probe only a small region of the 
inflation potential V{<p) and serve mainly to fix the value of the inflaton's dimensionless coupling constant {e.g. for 
V(4>) — A0 4 , we need A ~ 10 -14 ). Over-production of PBH's limits the value of Sh over many decades of "small" 
length scales (IMpc - 10~ 10 pc) and avoidance of such over-production, coupled with the COBE data, provides a 
powerful constraint on inflation |2fj|fl. 

Limits on black hole production are usually quantified by the parameter f3, defined to be the initial mass fraction 
of PBH's 

/3^4^« [ P(S)dS, (15) 
Ptot J § 

where (*) means we evaluate the ratio at the time of formation and the limits of integration are directly from (14). 
The probability distribution of density fluctuations, P(8), is usually assumed to be Gaussian, with the power spectrum 
giving us a = 5 rms = Sh at any particular scale. For the rest of this section we will assume P{S) is of Gaussian form 

"m-vk; «■*-£>• (16) 

and use the results as a gauge for any alternate distributions we find. 

There are two criteria for limiting the initial abundance of PBH's (see e.g. First, gravitational effects: the mass 
density of PBH's must not over-close the universe, SIpbh = {ppbh / ptot)\today < 1- And second, Hawking radiation 
from decaying PBH's must not disrupt any well-constrained physics {e.g. primordial nucleosynthesis). Since PBH's 
with M< 10 15 g will have decayed before today, gravitational effects for this mass range are non-existent. Similarly, 
larger mass PBH's have not had time to decay so Hawking radiation does not come into play. Thus, the physics of 
constraining PBH's is broken up into two mass regions: 

1. Gravitational constraints: M> 10 15 g 

2. Hawking radiation constraints: M< 10 15 g 
We discuss each in turn briefly. 

Limits from Gravitation: flpsn < 1 

We are interested in constraining the parameter 0, the initial density fraction of PBH's, using our limit on the current 
fraction: flpBH< !• We simply need to scale the density ratio with time. We consider PBH formation in the radiation 
dominated era, thus Qpbh ~ R until matter-radiation equality, after which the density fraction remains constant. If 
R* is the epoch of formation we have 



R* \ „ K ft 



1/2 



= Q PBH =l—jn PBH <lO- b \-j (17) 

where the inequality demands Qpbh< 1 and the time t* corresponds to R(t*) — R* . Since the black holes formed 
are roughly on the scale of the horizon at the time of collapse, we can write (15) in terms of mass scale using the 
following expression for horizon mass 

M H « 10 5 f - ] M (18) 



which yields 

0< 10" 



17 



M 
10 15 g 



1/2 



(19) 



The constraint becomes weaker for larger masses since higher mass PBH's form later and the density, ppbh, has a 
shorter time to grow relative to the total density. We now use (15), (16) to map the j3 values to limits on Sh- Note 
that since j3 is quite small, we depend very much on the tail of our assumed Gaussian distribution. For a 10 15 g PBH 
we have (3< 1(T 17 and the limit becomes 5 H < 0.04. Similarly, for M ~ 1M , S H < 0.06. 

Limits from Radiation 
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Hawking's theory of spontaneous black hole evaporation Q, J2lj predicts that black holes should emit particles at 
a characteristic temperature Tjj w 10 13 (Af/g) _1 GeV, which becomes more important with decreasing black hole 
mass. The lifetime of decay is r w 10~ 27 (M/g) 3 s p^], and is similarly dependent on the black hole mass, with the 
evaporation speeding up as the mass dwindles. PBH's of different initial masses evaporate at different times and their 
particle emission must not disrupt any well-understood physics associated with the time of their decay. We summarize 
a few constraints below. For more complete reviews see po|,^,^3| . 

• M ~ 10 14 - 10 15 g : Black holes in this mass range evaporate after recombination and must not produce a 7-ray 
density greater than the cosmic background [p4| . The constraint here is quite strong: PBH's could have at 
most 10~ 8 of the critical density today. We then have (3< 10~ 25 which corresponds to Sr < 0.03. We emphasize 
that since PBH production needs 5 ~ 1/3, formation in this mass range corresponds to a 10a effect assuming a 
Gaussian distribution of fluctuations. 

• M ~ 10 11 — 10 13 g: For these masses, evaporation occurs before recombination but the emitted radiation does 
not attain equilibrium due to the small value of the baryon to photon ratio. This effect will distort the background 
spectrum unless /3< 10~ 18 ( 1( jii g ) [ p5[ . For a typical mass, the result implies Sh< 0.04, and a Gaussian 8a 
fluctuation is required for PBH formation. 

• M > 10 10 g: This mass range marks the possibility of affecting big bang nucleosynthesis. For example, emitted 
photons from the evaporation must not photodissociate D. This limit gives pq |: (3< 10~ 21 ( j^o^ ) 1 ^ 2 an d 
corresponds to 8h< 0.03 for M ~ 10 10 g. PBH production again corresponds to a lOc Gaussian fluctuation. 

From the above discussions we see that PBH over-production provides limits on the initial spectrum of perturbations 
on mass scales small compared to the horizon mass ~ 10 55 g. We want to translate the above limits as a function of mass 
to limits as a function of co- moving length scale L. For lengths that cross inside the horizon before matter-radiation 
equality (L< 13/i _1 Mpc) the time of horizon crossing is 

t H ~ 3 x 10 8 (L/Mpc) 2 s (20) 
Together with (18) we have the following relation 

L x 2 



Mh ~ 30 UcV M& (21) 

with 1M corresponding to L ~ 0.2 pc and 10 15 g corresponding to L ~ 2 x 10 -10 pc. 

Since COBE fixes Sh ~ 3 X 10~ 5 on large scales, PBH constraints serve to limit the power spectrum from having a 
spike in power at some small length scale, or from having an excessively steep blue spectrum. Since Sh must be > 0.01 
in order for PBH formation to be important, we are concerned with factors of ~ 10 3 enhancement of the spectrum at 
small scales. We point out again that such a power spectrum is inherently associated with potentials that can make 
our Gaussian assumption invalid. This non-Gaussian tendency coupled with the fact that we are concerned with 
many-CT tails of the distribution (where even a small degree of skewness could drastically alter the results) suggests 
that any Gaussian analysis like the one above provides only a naive guess for the number of PBH's produced. We 
must understand something about the true fluctuation statistics before making any assumptions about the nature of 
the probability distribution far from the mean. 

III. STOCHASTIC INFLATION 



The stochastic analysis of inflation |27| provides an excellent method for obtaining statistics of scalar field fluctu- 
ations. In the stochastic approach, one divides the scalar field driving inflation, <p, into a long wavelength piece, 
and a short wavelength piece, 4> s 

4>{x,t) = <t>i(x,t) + <fi s (x,t) (22) 

where 

<t>l{x,t) = { <Pk6{tR{t)H -k)[a k tp k {t)e { - lk - x) +h.c] 



<j> s (x,t) = / d 3 k6(k - eR(t)H)[a k (p k (t)e(- ik ^ + h.c] (23) 



G 



and e < 1. The long wavelength piece contains only modes outside of the horizon: it is coarse-grained over the horizon 
scale and will eventually be interpreted as the "classical" inflaton. The short wavelength piece includes only modes 
inside the horizon, where quantum fluctuations are important. 

The goal is to determine how the sub-horizon quantum fluctuations affect the evolution of the coarse-grained field. 
We now briefly sketch the derivation of the stochastic slow-roll equation. (For a more complete review see [M and 
references therein.) First, insert the decomposition (22) into the Klein-Gordon equation, and write the result in the 
form 

j>l + 7^jV'(<j> l )=f(<j> s ;x,t), (24) 

where 

f(cf> s ;x,t) = ^ [4> s + 3H4> s -R- 2 (t) v 2 & 

As in the standard case (3), we have used the slow-roll approximation on 4n and have neglected the spatial derivative 
since this piece is homogeneous over the horizon scale. The <p s piece has its dynamics dominated by the spatial 
derivative and any effects of the potential are negligible. The high frequency modes are then approximated by those 
of a free massless field, ip k (t) , governed by the wave equation 

(^ + 3ff | + fc V 2 >^(i) = 0, (25) 

where we make the quasi-deSitter assumption: H « constant and R(t) ~ exp(iit). The solution to the free field 
equation is known 

<p%(t) = -J= (77 - i/k) exp(-ife77) (26) 
17 : 



2k 
dt 



R(t) H ' 

so we have an explicit expression for <p s by letting <^fc(t) = <fi k (t) in (23). Inserting the expression for <f> 3 into (25) we 
have 



/(M) = g 



o 2 _„a 



3Hj- + e - 2Ht \/ 2 



dt 2 dt 
ieR(t)H 



d 3 k9{k - eR(t)H)[a k ip* k (t)e ( - tk ^ + h.c 



V2(27rfc) 3 / 2 



d A k5{k - eR(t)H) [a k e~ tk - x - h.c] (27) 



and equation (24) is now complete. The field 4>i is smoothed over the horizon scale and we are concerned with the 
dynamics at a single spatial point. Let <pi(x,t) — > 4>i(t) — > <j>(t) and interpret this piece as a classical field which is 
acted on by the stochastic force f(x,t) — > f(t). The field <fi is affected not only by V , as in the standard slow-roll 
approach, but also by the flow of initially small-scale quantum fluctuations across the horizon. If we calculate the 
expectation value and two-point function of f(t) and interpret them classically we obtain 

(/(*)) - 0, 

{f(t)f(t')) = ^5(t-t>), (28) 

and (24) becomes a Langevin equation for (j>. Regrouping terms and slightly changing notation, the final stochastic 
slow-roll equation is 

1 #3/2 

V'(cj>) + ?—g(t) (29) 



ZH{4>) ' 2tt 
where g(t) has the properties of Gaussian white noise 



(<?(*)> = 0, 

(g(t)g(t'))=6(t-t'). (30) 
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We have standard slow-roll evolution with an additional stochastic term which represents the effect of quantum 
fluctuations on the dynamics. 

We will use the stochastic equation (29) to determine the probability distribution of fluctuations, P((j>, t), for models 
of interest, and then measure the degree to which non-Gaussian statistics are important. The Fokker-Planck equation 
associated with P(4>, t) and (29) is well known [M 



d t P( 



3H(cj>) 



P(0,t) 



#3/2 



(0)d (p 3 / 2 (0)P(<M) 



(31) 



with the Stratonovich interpretation of the noise. 

Several authors have employed the stochastic Langevin equation and corresponding Fokker-Planck equation to 
explore the importance of non-Gaussian fluctuations for large-scale structure po| , p8|JrI| . We use a similar approach 
but on the much smaller scales relevant for PBH formation. We also choose to work with the more intuitive Langevin 
equation (29) when determining our solutions rather than the more complicated Fokker-Planck mathematics of (31). 
Next we present several examples which illustrate our method of solution and provide a feel for how non-Gaussian 
statistics arise from non-linear inflaton dynamics. 

• de Sitter: The Gaussian Case 

For our template example we examine the case of constant vacuum energy, where there is no possibility of mode- 
mode coupling and the statistics should be exactly Gaussian. To show this point clearly, start with equation (29). 
We have H((f>) = H = constant and V'(4>) — 0, so the evolution is simple diffusion 



P 3 / 2 
2tt 



g(t). 



The normal statistical behavior becomes apparent after integrating 

#3/2 rt 



<P(t) = <j} 



2n 



g(t')dt'. 



(32) 



(33) 



Since g(t')dt' is Gaussian distributed as a noise source, the resulting probability distribution for <j> must be Gaussian 
as well. From (33) we have 



W) - 



= o, 



4ir 2 



and the probability distribution 



P{4>, t) cx exp 



m) - ci> f 

2a 2 



(34) 



(35) 



a 2 ^(\^t)-4> \ 2 ) = : —{t-t ) 



IP 

4-7T 

Over a Hubble time (t — t a = H^ 1 ), u = H/2tt: the standard result from quantum field theory. 

Non-Gaussian distributions arise when the evolution is non-linear, or when the inflation is not strictly de Sitter. 
One gains insight into the statistics of more complicated examples by comparing them to the simple de Sitter case. 
Our goal when finding distributions for our examples will be to put the stochastic equation into the form of (32) so 
that we may simply write down the answer using (35). We illustrate this technique in the examples that follow, and 
use the notation: 



(j)/m p f, t-*t/m p i 



in order to keep the variables dimensionless. 
• Non-linear Diffusion: Driftless d> 4 
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Mode-mode coupling arises when either the V "drift" term or the diffusion term in (29) is non-linear in <j>. In order 
to develop some intuition for how each of these effects generates non-Gaussian statistics, consider first an example 
with non- linear diffusion when the drift component is ignored (V — > 0). For concreteness let 



V 



X 



4 4 



(36) 



with H 



'2ttA „,2 
3 



x 2 . If wc ignore the drift term in (29) our stochastic equation is 

x = Cx 3 g{t), (37) 
Z(x), such that Z(x) will satisfy the simple 



where C = A 3 / 4 [547r]~ 1/4 . The trick now is to change variables, x 
diffusion equation (32). We want: 

Z = = \px 3 g{£)\ — [constant]g(t), 
where we have used (37) in the second step. Clearly if we let 

then we are left with an equation of constant diffusion 

Z = Cg(t). 

Now wc write down the answer in correspondence with (35) 

(Z 



P(Z, t) (x exp 



Z Q f 



2CH 



Finally, changing back to x yields the probability distribution of interest: 



P(x,t) 



dZ 



dx 



P(Z(x),t) (xx 3 exp 



2\2 



(s-a-O 



8C 2 t 



(38) 



(39) 



(40) 



(41) 



(42) 



The statistics here are clearly non-Gaussian, but is the deviation significant? The answer depends on the value 
of C 2 t which is, in an approximate sense, the effective width of the distribution, cr e ff- Expression (42) is plotted in 
Fig. 1 for C 2 t — 1 and xo = 1 with arbitrary normalization. We see that P(x) has a long tail of large fluctuations, 
and clearly is far from Gaussian. For small values of cr e ff, the deviation of x from x will not be significant and the 
distribution will recover a Gaussian shape. We will discuss this behavior in detail following the next example. 

• Non-linear Diffusion: Full 4 Theory 

Now let us go one step further and include the diffusion term for the potential (36). We have V' = Xm 3 ^ 3 and the 
stochastic equation (29) becomes 



x = -C\x + C 2 x 3 g(t). 



(43) 



where C\ = and C 2 = A 3 / 4 [54tt]- 1 /4_ xhe drift term in this 

case is linear and does not lead to mode-mode 
coupling. The probability distribution of fluctuations will then look very much like that of the driftless case (42), 
except that the peak will shift to smaller x values as the field rolls down the hill, and the spread in the distribution 
with time will no longer be linear. 

As before, we obtain the solution through a change in variables: x(t) — > Z(x,t). In this case we would like to first 
remove the drift term by effectively going to a frame moving with the mean "classical" velocity (thus, the explicit 
time dependence in Z). By the classical velocity we mean the solution to (43) when the quantum fluctuating is term 
left out: 



x = —C\x, 



(44) 



which is simply 
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xa(t) = x exp(-Cxt). 

We would like our new variable to be constant as long as x(t) follows this classical path, so we let 

Z(x,t) — xexp(C\t) 



x c i(t) 



We have chosen Z{x = x c i{t), t) — x Q . Now Z has time derivative 

Z = xexp(Cit) + Cixexp(Cit), 
and using (43) for x we have a stochastic equation for Z 

Z = C 2 x 3 exp(Cxt)g(t) = C 2 Z 3 exp(-2Cit)g(t). 



(45) 



(46) 



(47) 



(48) 



Now our equation is much like (37) except for the extra time-dependent factor. We can, however, scale this factor 
out using a change in time variable: t — > r(t). If we remember that g 2 (t) ~ 6(t) then we have 



37 

— = C 2 Z 3 exp(-2C 1 ) 



dt 

rf7 



-1/2 



9(r). 



In order to match (37) we need dr/dt — exp(— AC it), and requiring r(t = 0) = gives 

r(t) = [1 " exp(-4C 1 t)] = ±- [l - (x cl (t)/x ) 4 ] . 

We now havefj] 

Z = C 2 Z 3 9 (t) 

just as in (37), and we simply use the solution (42) to give us the probability distribution for Z 

(z-*-z-*r 



8Cfr 



P{Z, t) cx Z' 3 exp 

L 0{ ~"t 

We transform back to x and t and obtain the standard result for </> 4 theory 

(.T- 2 -:r c/ (£r 2 ) 2 



(49) 



(50) 



(51) 



(52) 



P(x, t) cx x exp 



(g) [{x /x cl (t)f 



-1 



(53) 



So, as expected, the distribution has the same non-Gaussian form as did (42). In truth, however, most practical cases 
[i.e. those concerning large scale structure) find this distribution quite well approximated by a Gaussian due to its 
extremely small standard deviation, a oc 

To put in some numbers, let us assume we are interested in the statistics of fluctuations associated with large- 
scales, say the modes that cross outside of the horizon during the epoch x ~ 4. We need to evaluate (53) at 
x c i(t) = x en d ~ 0.4 in order to find the distribution of interest. Since for this model A ~ 10 -14 , let us assume that the 
standard size of a fluctuation about x en d in the final distribution is quite small: x = x en d + e (where e/x is assumed 
Cl). In (53) we have 



X end 2e/x end 



and 



P(e) cx exp( 



2a 



(54) 



(55) 



From now on, the over-dot will signify the derivative with respect to the argument of the stochastic function g. 
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which is Gaussian as expected with an effective standard deviation 



e 2 // = - J %f [(xo/x end r " 1] = 12 [(Waw) 4 - !] ^ (56) 



Plugging in the numbers we have cr e // ~ 2 x 10~ 7 . So our approximation e/x en d <C 1 is clearly self-consistent, and 
the Gaussian assumption is a very good one in this typical case. 

The small value of <x e / / is responsible for the near Gaussian statistics, but the more the mean size of fluctuations 
increases (<J e // T)> the more (54) fails and the Gaussian approximation (55) breaks down. In this way we can 
understand why a flat potential region gives rise to non-Gaussianity. Let us use C\ as a crude measure of the 
importance of the drift term (<~ V') in the dynamics (43). Notice that as the drift term gets smaller (C\ j), the 
effective width of the distribution (56) grows, and with it the likelihood of non-Gaussian statistics. Notice also that 
as we consider fluctuations farther and farther from the mean, approximation (54) will get worse and worse. The 
number of standard deviations one must be from the mean in order for non-Gaussian statistics to be important is 
very sensitive to the value of A. However, we can at least see the qualitative effect that large fluctuations have on 
highlighting any intrinsic non-Gaussianity. 

• Non-linear Drift 

Finally we examine the case of non-linear drift as a source of non-Gaussian behavior. Consider the potential 

V = Xrripi(l + ax 3 ) (57) 

over a region where |x| <C (1/a) 1 / 3 . The potential is roughly constant over this region (V w Am 4 ;) so the only 
non-linearity comes from V = 3\m 3 l ax 2 . 

The stochastic equation in dimensionless variables is then 

i« -Cix 2 + C 2 g{t), (58) 

where C x = a^3A/87r, and C 2 = A 3 / 4 (32/27?r) 1/4 . Since the diffusion term is now constant, any non-Gaussian 
statistics will arise due to the x 2 drift term in (58). As before, we solve for the probability distribution by making 
a change in variables, x{t) — ► Z(x,t). We want to remove the drift term by changing to a frame moving with the 
classical velocity and canceling out the drift in (58). The classical solution to the standard slow- roll equation of 
motion 

x = -Cix 2 (59) 

is trivial 

ar c j(*) = (Cit + a:- 1 )- 1 ) (60) 

where x a = x c i(t = 0). Let us pick x a = —1 for concreteness. Our new variable, Z(x,t), should be fixed as long as 
x(t) travels along its classical path: -^Z(x c i,t) — 0. We achieve this goal by letting 

Z(x,t) = {x- 1 - at)' 1 , (61) 

where we have picked Z(x c i(t), t) — x D . Differentiating with respect to t and using (58) we obtain a driftless stochastic 
equation for our new variable 

Z = %C 2 g(t). (62) 
x £ 

At this point let us simplify the problem by approximating Z and x with their mean classical values: x — > x c i (t) and 
Z -» Z(x c i(t),t) =x D = -1. We have 

Z = (at + lfag(t) (63) 

which is beginning to look a lot like simple diffusion. Clearly Z(x,t) will be Gaussian distributed, but <r(t) will be 
more complicated than in (35) due to the explicit i-dependence in the diffusion coefficient. To find the correct form 
of a(t) let us change time variables t — ► r such that the time dependent diffusion will be scaled away. Remembering 
that g 2 (t) ~ S(t) we have 
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dZ 



So we need dr/dt — (C\t — l) 4 , or 



[Cii - if C 2 



1 



~dr~ 
~dt 



-1/2 



5(t). 



r = ^[l + (^-l) 5 ] 



Then our final stochastic equation for Z(t) mirrors (32) 

= C 2 g(r) 

and corresponds to the probability distribution (35): 

{Z + lf 



dZ_ 
dr 



P(Z, t) ex exp 



2ct(t) 2 



■o{rf=ClT. 



We obtain the distribution of interest, P(x,t), by simply changing back to the appropriate variables Z — > x, r 
We have 



P(x,t) cx 



8Z 




dx 


exp 1 



2cq T (t) 



(64) 

(65) 

(66) 

(67) 
-> i. 

(68) 



Note that the non-Gaussian behavior is evident not only in the non-linear a;— dependence of Z, but also in the non- 
linear time evolution of r(t). If we evaluate the distribution at a particular time, say that corresponding to x c i = —2, 
we will have a peak in probability near this value. From (44) we need t = (2Ci) _1 , which corresponds to r « (5Ci) _1 
from (65), so (68) becomes 



P{x) oc — J exp 



([x- 1 -0.5]- 1 + 1 
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2ct 2 



(69) 



where a 2 sa C 2 2 /5Ci. 

We plot (69) for <r 2 = 0.02 in Fig. 2. The non-Gaussian behavior is clear. There is a long tail in this distribution for 
small fluctuations and a sharp cutoff for large fluctuations — in opposition to the previous example. The non-linear 
drift tends to encourage fluctuations to fall down the hill (since V" < the slope gets steeper for smaller values of x) 
resulting in an under-abundance of large fluctuations compared to the Gaussian. 



IV. TOY MODELS 



In this section we present three example potentials which meet the small-scale power requirement for PBH produc- 
tion (Sh> 0.01) as discussed in Sec. II. For each example we apply the methods developed in Sec. Ill to determine 
the probability distribution of fluctuations and compare this distribution to the usual Gaussian assumption. Each 
example is realistic in the sense that it forces the power on large scales to be small (6h ~ 3 x 10 -5 ) while giving 
us the sufficient small-scale power needed. The additional constraint on the inflaton potential is associated with the 
over-production of gravitons and limits the scale of inflation: H/m p i< 10 -5 during the epoch associated with the 
horizon scale today (</>(10 4 Mpc)) The first two models below are actually on the hairy edge of this constraint, 
but recall that we are not trying to construct "true" models here, we are simply interested in the non-Gaussianity 
associated with PBH production. Since the non-Gaussian statistics in these models derives from the shape of the 
potential, slightly altering the scale of inflation will not drastically affect the result. The gravitational wave constraint 
on our third toy model is a different story and we will discuss the reason when we come to it. 

We emphasize that these toy potentials are meant only to illustrate the importance of non-Gaussianity in PBH 
producing models. The correct distribution for any "realistic" model will depend on the nature of the potential. The 
idea is that the following examples will be tools of pedagogy, enabling the reader to easily perform calculations of 
PBH abundances without relying on an incorrect Gaussian assumption. Before we begin, let us restate 

x = 4>/m pU t—tt/rripi, 
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and introduce some notation 

V = XmfaV, H = VXm p i/H. 

The symbols V and H are dimensionless quantities of order unity which will be useful in keeping track of small 
parameters. 

A. Plateau Potential 

Consider the potential, shown in Fig. 3, which in some region of interest (x ~ 0) has a flat "plateau" feature ^ 

T %. _ f 1 + arctan(x) x > 

^(l + (4 x lO 33 ^ 21 x < 0. [W > 

We have designed this somewhat outrageous potential to produce a spike in small-scale power using the methods 
outlined in Sec. I. It suffices to regard expression (70) as only a region of potential, so we arbitrarily begin following 
the roll-down of x at the value corresponding to the scale L ss 10 pc, or x st ~ 0.15. We choose to normalize our 
potential region at this hand-picked starting point, and fix A with <5ff (10 pc) = 3 x 10 , which is an arbitrary but 
reasonable choice. Expression (4) then requires A s» 6 x 10~ 10 . When the field first reaches the plateau region, it 
is moving too fast to obey the slow-roll conditions, and we are forced to obtain the form of 8h by solving (1),(2) 
numerically. We also use the more exact behavior 8h oc H 2 /x. The plateau region is so flat, however, that the friction 
of the expansion slows the field down very quickly. The magnitude of Sh grows as x slows down and continues to grow 
until x reaches its lowest value, corresponding to a peak in power, and also corresponding to the beginning of a second 
slow- roll regime. (Note that this potential provides two periods of inflation with only one inflaton.) Numerically, we 
find that the peak in power occurs when = —1.23 x 10~ 2 . The height of the peak is then 

5 H (x») ~ VA^-O*) ~ 0.05. (71) 

We show the density perturbation spectrum at horizon crossing associated with this region in Fig. 4, and see that 
the peak in power corresponds to the creation of ~ IMq black holes. 

The modes which will eventually be responsible for PBH production pass outside of the horizon at the epoch x*. 
Since we are only interested in calculating the non-Gaussian behavior associated with these modes, we evolve the 
stochastic equation from x* — > x enc i to determine the relevant probability distribution. For this model, the end of 
inflation occurs when the first slow-roll condition breaks down: 



v" 



2iirV 



1, and corresponds to x en d — —1.55 x 10~ 



The stochastic equation (29) in the new notation is 

-VXV' A 3 / 4 ff 3 / 2 

x = = 1 

3H 2tt 



»(*)• (72) 



Note the relative sizes of the drift and diffusion terms in the evolution of x. In most cases, the higher power of A in 
the diffusion term means that quantum fluctuations play a very minor role in the dynamics of x, resulting in nearly 
Gaussian statistics. But in this example, the drift term (V) is very tiny over the region of interest, and the two terms 
become of more equal weight, increasing the likelihood of non-Gaussian statistics. 

Over the range of interest, |x| ~ 10 -2 , so to a high degree of accuracy, V ~ 1, and we have a case of constant 
diffusion. On the other hand, the drift force is small and very non-linear: V' — 8.4 x 10 34 x 20 . This situation is very 
much like the last example presented in Sec. Ill and the method of solution will be similar. Let us re-scale x by x* to 
keep our variables of order unity: x = x/\x*\. The path of the mean evolution is then from x„ = —1 to x en d — —1.26. 
Using A = 6x 10~ 10 in (72) we have our stochastic equation for x 

2 = -(1.2 x 10" 7 )i 20 + (7.8 x 10- 6 )g(t). (73) 



This potential with a single break is similar to the double-break potential proposed by Ivanov, Naselsky, and Novikov 
for making PBHs. 
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We see that the small value of V' gives the quantum diffusion term a more equal weight in the dynamics. Let us do 
one more bit of cleaning up by rescaling the time variable t — > t(1.2 x 1CP 7 ) (see (49)) which gives 



x = -x 20 



ag(t') 



(74) 



a = 2.3 x 10" 



We now mirror the method of solution presented in Sec. III. 

The first step is to go to a frame moving with the classical velocity. The classical slow-roll equation is x = i 20 and 
using i = x* = — 1, the classical path is 

x cl (t) = -[l-19t}- 1 / 19 . (75) 
Our variable, Z(x, £), will remain constant along x = x c i 

Z(x,t)^[19t^x~ 19 }^/ 19 (76) 
where we have chosen Z(x c i(t),t) = —x Q = 1. The stochastic equation for Z is then diffusion-only by design 



20 

Z = - [ | ) ag(t), 



(77) 



and approximating x — > x c i(t), Z — > — x we have 

Z~ (l-19i) 20/19 ag(i). (78) 
Now, with one more change of variables we can scale away the time dependence multiplying g(t). Let 

r = 1 (l - (1 - 19£) 59 / 19 ) , (79) 
so that we are left with an equation exactly like (32) 

Z = ag(r). (80) 

The probability distribution, P{Z, r), follows from (35), and we change variables back again to obtain the distribution 
of interest: 



P(£, t) oc exp 



([19i - ^- 19 ]- 1 / 19 - I)' 
2a 2 r(t) 



(81) 



We want to evaluate the distribution at t ene [, the time corresponding to x c i = —1.26. From (75) we find 19t enc i — 0.9876 
and from (79) we have T(t en d) ~ 1/59. Plugging in these values along with a = 2.3 x 10 -2 , we have the final distribution 
of fluctuations: 



P{x) oc -^tt exp 



(i 



K 19i e 



-1/19 



1 



(82) 



where K — 5.6 x 10 4 . The distribution is plotted in Fig. 5 along with the results of a numerical simulation for this 
same potential. For our simulation, we started with the stochastic equation (72) and made no approximations. We 
used the Box-Muller method |]33|l to transform uniform deviates into random Gaussian deviates to mimic the stochastic 
force. The numerical results consist of 4 x 10 4 individual runs of the Langevin equation, and we have normalized the 
height of our calculated distribution to fit this number. We see that the calculated distribution (82) agrees well with 
the numerical results. 

The distribution of fluctuations is clearly non-Gaussian. In Fig. 6 we have plotted our analytic distribution (82) 
along with two Gaussians that one may wish to compare it to: one with the same mean and standard deviation as 
our distribution, and one with the same height and width at half maximum as the peak in our distribution. We see 
that the true distribution is under-producing large fluctuations compared with either of the Gaussians. Now, for PBH 
production under the Gaussian assumption, we are concerned with fluctuations on the tail of the distribution ~ 6<r. 
The calculated distribution differs so drastically from the Gaussian assumption at this distance from the mean, that 
we can only compare them on a log scale. Fig. 7 shows distribution (82) along with the two Gaussian comparisons 
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in units of the standard deviation from the mean. As the more conservative choice, we use the a associated with the 
Gaussian fitted to the peak of (82). Observe that the Gaussian assumption would vastly over-produce PBHs, with 
an error of order ~ 10 150 at 6<t! This model was designed to give us a large number of PBH's using the Gaussian 
assumption, but upon examining the distribution more closely, we see that the actual production is practically zero. 
We, as inflation designers, are forced to make an even higher spike in small-scale power if we want PBH production. 
Note also that as the spike becomes higher, the drift term becomes smaller, and the distribution will tend to skew 
even more towards small fluctuations. 



B. Wiggle Potential 



Our second potential which produces significant small-scale power is one with a wiggle in the path of the inflaton 

V = 1 + [136.717]x 3 - 0.05x (83) 

as shown in Fig. 8. We start our evolution at x(L = 10 4 Mpc) = 2.6, and since slow-roll is obviously invalid over 
the dip, we integrate (1),(2) to obtain the spectrum of 8h{L) shown in Fig. 9. Normalizing at x = 2.6 to COBE gives 
A~5x 10~ 14 . We have adjusted the wiggle to make the field slow down dramatically just after the top of the bump 
(xtop = —1.0411 x 10~ 2 ) which produces a large spike in power corresponding to the ~ 10 28 g mass scale. We find 
the slow-point to be x* = —1.10448 x 10~ 2 , which is also the beginning of a second slow-roll epoch of inflation. The 
height of the peak is Sh ~ */\^j-(x*) ~ 0.01, the order of magnitude we need for PBH production. 

After the field slows down at a;*, inflation continues until x en d ~ —0.1. This is the path of interest for estimating 
the distribution of fluctuations. Let us again set our notation before writing down the Langevin equation. For the 
analytic calculation we use V ~ 1 over the range of interest^] and also scale x by x — x/xt op , such that x% op = — 1 
corresponds to the top of the wiggle (V'(x = —1) = 0). Using this notation, the starting point is x* — —1.00033 and 
we want to follow the evolution until x enc i ~ —9. The potential in terms of £ is 

v re i, v' = om[x 2 - 1] 

and the stochastic equation for x follows 

x = -(1.2 x 10~ 7 )(i 2 - 1) + (7.5 x 10~ 9 )g(t). (84) 

As before, we clean things up with a new time scale t — > £(1.2 x 10 -7 ), which gives us 

£ = -(x 2 - 1) + ag(t), a = 2.2xl0" 5 . (85) 

We want to use the same trick and factor out the drift term by changing to a variable which is constant along the 
classical path. We solve £ = — (i 2 — 1) to obtain the classical solution with the initial value £ c i (t = 0) = £ tf — —1.00033 
and obtain 

£ cl (t) = coth(i - 4.35). (86) 
Following the methods outlined in Sec. Ill, the new variable, Z(x,t), is 

Z(£, t) = coth[coth _1 (i) - t] (87) 

where 

Z(x cl (t),t) = x*. (88) 

The stochastic equation for Z is then 

Z = a%f^g{t). (89) 
x z — 1 



6 This approximation is very good at x* but is off by ~ 10% at the end of inflation. However, since the field spends most of 
its evolution time near x», this approximation should be fine. 



15 



If we use the approximation x — > x c i (t), Z — + we have 



coth^i - 4.35] - 1 



= sinh 2 [i- 4.35] bg(t), 

where b = 1.47 x 10 -8 . Now we want to change to a new time coordinate, r(t), which will leave us with a simple 
diffusion equation. The requirement is 

^ =smh 4 [t-4.35l, (90) 
at 

or 

r(t) =C+\[t- 4.35] - - sinh[2(i - 4.35)] + ^- sinh[4(i - 4.35)], (91) 
8 4 32 

where C — 5.6 x 10 5 demands that r(t = 0) = 0. So the variable Z(x,t) obeys simple diffusion by design 

Z = bg(r) (92) 

and has a Gaussian probability distribution with mean Z = 2* and a 2 — b 2 r. Now, as discussed earlier, we simply 
perform the reverse transform Z — > x, r — > t to obtain the distribution of inflaton fluctuations: 



« ( ~^tt 1 exp 



(coth(coth - t) - i*) 2 ] 



2b 2 r(t) 



We want to evaluate the above expression at the time t en ^ — 4.24, when the classical path reaches x en d = —9. 
Equation (91) gives us r(t en d) ~ C and plugging in all of these values we have the final probability distribution 

P(x) cx - j exp -K (coth(coth _1 (i) - 4.24) + 1.00033) 2 , (94) 

with K — 4.1 x 10 9 . As in the previous example, we have simulated this distribution numerically using no approxi- 
mations, and we plot the two together in Fig. 10. We have normalized (94) to the simulation height. We see again 
that the derivation does quite well, and the distribution has a deficit of large fluctuations relative to a Gaussian. The 
distribution is skewed negative again since the drift term ~ V and V" < 0. Negative fluctuations tend to fall down 
the hill more quickly than positive ones. There is no need to do another explicit comparison to a Gaussian distribution 
since if the distribution is clearly non-Gaussian to the eye, then the high-er tail will be exponentially worse. 

C. Cliff Potential 

For our last example we present a potential region which flattens out to achieve large power in the form of a blue 
spectrum, and then has a cliff-like feature where the slope abruptly becomes much steeper. Consider the region of 
potential, shown in Figure 11, given by 

4 / cos- 2 [1.5a:] x > 0.00004 , , 

V ~ p 1 \ {2.7)- 4 [x + 2.7] 4 x < 0.00004. 

For x > 0.00004 we use a form suggested by Hodges and Blumenthal jl(| which gives us Sh ~ (L/L )~ n , n ~ 0.18. 
The highest amplitude in power, corresponding to the flattest region of the potential and PBH formation, occurs at 
= 0.00004, after which the slope increases abruptly as a " 4 " form ends inflation. Since we are interested in 
fluctuation statistics which depend only on the nature of the potential from x* — > x enc i, the "small V"' argument 
does not apply to encourage non-Gaussian statistics in this case. Using A~3x 10 -11 we obtain the power spectrum 
shown in Fig. 12. The highest point in density is 5h ~ 0.04 and corresponds to ~ 10 20 g PBH's. 

As we mentioned before, the region of interest for calculating the probability distribution has quartic form. Recall 
that we have already calculated the distribution for the quartic case in Sec. Ill (53) so our job is nearly done. In 
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order to make things look more familiar, define y = x + 2.7, A' = 4A/ (2.7) 4 ~ 2. x 10 12 , such that we have the more 
standard form 

V(y) = jy\ (96) 

and we are interested in the path from y ss 2.7 to y en d ~ 0.4. Now, letting A — > A' and x — > y in (53) we have the 
fluctuation distribution we need 



P(y) oc y 3 exp 



dr 2 - u~ 2 ) 2 ~\ 

■' ')■/ ' (97) 



2a 2 



0-2 = -j[(y*/yend) 4 - l]. 

From the discussion following (53), we know that P(y) is skewed positive, but we also know that the amount of 
skewness depends on the value of a 2 . In this case, non-Gaussianity is negligible due to the extremely tiny value of 
A'. In order to see the lack of skewness in a quantitative way, let us borrow a measure proposed by Yi, Vishniac, and 
Mineshige jll). This simple estimate for skewness is the ratio 

R = P ^ d + Na ^f\ „ exp \N*Vwr±(yi - yL d ) 1/2 } « exp [^(10~ 5 )] (98) 

P(y e nd- Na e ff) L J 

where cr e ff, the number of effective standard deviations from the mean, was given in (56), and N measures the size 
of the fluctuation. We see clearly that the small value of A' tends to force R — > 1 and prevent skewness. But a large 
value of N (e.g. a large fluctuation) will compete with this effect and give a skewing effect. For our potential (94) we 
have chosen the largest value of A — > A' possible to be consistent with the gravitational wave constraint H ~ 10 -5 . 
However, even with the large fluctuation we need for PBH formation, N ~ 10, expression (97) implies that the skew 
is only 1%. Suppose, for example, that we did ignore the gravity wave constraint Q and push the value of A' up to 
~ 10~ 9 . The ratio (97) would then give us a ~ 20% effect at N — 10. Thus even for potentials with "normal" shapes, 
non-Gaussianity may be important due to the rare large fluctuations associated with PBH production. 

V. CONCLUSION 

We have argued that the very conditions associated with PBH formation: 

1. Large spikes in power, usually associated with flat potential regions 

2. Rare, many- a fluctuations 

are exactly the same conditions which can produce significant non-Gaussian fluctuations. Flat potential regions 
promote the importance of quantum fluctuations in the inflaton dynamics and encourage the mode-mode coupling 
responsible for non-Gaussian statistics. On top of this effect, many-cr fluctuations push us out to the tail of the 
probability distribution, where any intrinsic skewness will be amplified. We have quantified this intuition with several 
toy models that produce the small-scale power associated with PBH production, and have used the the stochastic slow- 
roll equation to obtain the fluctuation distributions. Our examples clearly illustrate that the Gaussian assumption can 
lead to large errors in the calculated number density of PBH's, and that the nature of the non-Gaussian distribution 
is extremely model dependent. 

Specifically, for models with spikes in small scale power, the fluctuation distributions were skewed towards small 
fluctuations, under-producing PBH's by many orders of magnitude relative to the Gaussian assumption. The negative 
skewing in these examples came from mode- mode coupling due mainly to non-linear drift, which encouraged negative 
fluctuations to fall down the hill faster than positive ones (V" < 0) 0. Because the fluctuation statistics depend on the 
path of the inflaton after it passes the flattest region of the potential, we were even able to construct an example where 
the Gaussian assumption holds, simply by forcing a dramatic increase in V' (a "cliff") just after the peak in power. 



7 Note that we can easily adjust our potential region to obtain the appropriate small scale power with a new value of A. 
8 A potential region which tends to "cup" the inflaton (V" > 0) after the flat region may produce positive skewing. 
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One may regard such a cliff region as unnatural, but it does illustrate the model-dependent nature of distribution 
shapes. 

These results have several important implications and we discuss each briefly. First, the standard approach for 
limiting the initial fraction of PBH's, /3, and limiting the spectrum of initial density perturbations (see Sec. II) must 
be reconsidered. Because of the model-dependent nature of the distributions, it is not possible, as in the standard 
practice [Q, to use PBH over-production to limit generic inflationary power spectra. Correspondingly, it is not 
possible to determine the number of PBH's produced from the power spectrum alone. We can use our toy models to 
exemplify what might happen in various cases. For models with spikes in power, like our plateau and wiggle examples, 
fluctuation distributions will probably skew towards small fluctuations and under-produce large fluctuations relative 
to the Gaussian case. Many-er fluctuations then will be much less likely, and PBH production in these models will 
require an even higher spike in power. Then magnitude limits on spikes in power to prevent PBH over-production 
will be less stringent than the limits obtained using Gaussian statistics. We also point out that the formation of 
any PBH's at all without over-producing them will require an even more drastic fine-tuning than in the Gaussian 
case. The fine-tuning must be more precise because the large fluctuation tails in these models are much steeper than 
a Gaussian tail. So small changes in fluctuation size will cause a larger change in the probability of having such 
fluctuations. 

Our results also affect the use of PBH over-production to rule out or limit the parameter space of specific inflationary 
models. For example, authors employ PBH over-production to constrain the slope of blue perturbation spectra from 
inflation (see J23|). All such examinations use the Gaussian assumption, and limit models without regard to the shape 
of the potential after the power reaches its maximum height. But the amount of PBH production, and hence constraints 
on 5h and the slope of the spectrum, depend crucially on the nature of the fluctuation distribution and hence on the 
shape of the potential between the region of PBH formation and the end of inflation. Our "cliff" potential toy model 
actually recovers the Gaussian approximation, but only because the slope of the potential increases dramatically just 
after the region of high power associated with PBH production. However, with a more rounded potential shape 
after the peak instead of a cliff, non-Gaussian fluctuations will be much more likely. PBH over-production is also 
important for constraining parameter space in hybrid inflationary models due to associated spikes in small-scale 
power J5|,|| . Again, our toy models indicate that spikes in power are associated with PBH under-production relative 
to the Gaussian case. However, without further investigation of multiple-field PBH production models, the validity 
of our intuition from single-field examples remains unclear []. 

Authors occasionally investigate the possibility of PBH formation associated with a soft equation of state during 
ap = "dust like" epoch p4| , p3| . Dust era PBH formation would be important if the universe underwent an early 
p = stage (before the usual matter-dominated epoch) . Conceivably, such a dust stage could occur due to some as yet 
unknown physics (e.g. possibly due to coherent inflaton oscillation during reheating j35|). For dust era PBH formation, 
the major criterion for PBH collapse is that the initial perturbations be non-rotating and spherically symmetric p4] , p3"f , 
and the probability for spherical geometry typically increases with the amplitude of a fluctuation. Again the associated 
power spectra must have excess power on small scales to achieve substantial PBH production, and the fluctuation 
distributions will typically be non-Gaussian. But the tails of probability distributions are not so important to the 
analysis since rare large-amplitude fluctuations are not all that is important here, but also spherical ones. Thus 
the value of (3 will depend mainly on (the rms) 8h, with only a weak dependence on the shape of the distribution. 
Compared to PBH formation during the radiation dominated era (where, as we have shown, the distribution is very 
important and non-Gaussianity can shift (3 by many orders of magnitude), the effect of non-Gaussianity on dust era 
formation should be much smaller. 

Finally, we discuss our results in the context of PBH formation associated with the QCD phase transition. 
Jedamzik M points out that PBH formation due to a first order QCD phase transition would be roughly consis- 
tent with ~ 0.5Mq MACHO's, since this is roughly the total mass-energy inside the horizon during the QCD epoch. 
However, even if the QCD transition is first order, the minimum 5h for PBH formation at this mass scale will be 
lessened by at most a factor of order unity, for the following reason. During the epoch of the first order phase transi- 
tion, the effective velocity of sound drops to zero, and with it the Jeans mass. Thus density perturbations which cross 
inside the horizon at the beginning of the epoch can begin to grow. But the duration of the phase transition is quite 
short. Fluctuation amplitudes on this scale will be enhanced slightly relative to the standard case before the universe 
again achieves a hard equation of state. If the fluctuation amplitude at this time is as large as S ~ 1/3 then PBH 
formation can occur. So the minimum value of the initial Sh at horizon crossing is reduced, but only by a factor of 
order unity (see |56f for a more complete discussion of perturbation growth during this epoch) . This slight increase 



9 An analytical examination of multiple field models may prove to be too difficult, but numerical simulation is always possible. 
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in the PBH mass function would be important if PBH production were marginally possible over some range of mass 
scales. PBH formation would then be enhanced for masses around ~ 0.5Mq, which could perhaps explain why this 
mass range of MACHO 's is observed. But in order to achieve PBH formation at the QCD mass scale, we still must 
have large-amplitude fluctuations at this wavelength, so our arguments for non-Gaussianity still apply. Moreover, any 
non-Gaussianity in the fluctuation distribution should be important, since PBH formation is again associated with 
rare fluctuations, and therefore quite dependent on the shape of the distribution far from the mean. 
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FIG. 1. The unnormalized probability distribution (42) for the driftless </> 4 example, where x = <f>/m p i. We use the the 
value x = 1 arbitrarily and pick C 2 t — 1, which is unnaturally large, in order to demonstrate the intrinsic non-Gaussianity 
resulting from non-linear diffusion. 



FIG. 2. The unnormalized probability distribution (69) for the non-linear drift example. We have used a 2 = 0.02. Note that 
the non- linear drift has caused negative skewing, in contrast to Fig. 1, due to the tendency for negative fluctuations to fall 
down the hill more quickly than positive ones. 



FIG. 3. The plateau potential (70): V(x > 0) = 1 + arctan(a;), V(x < 0) = 1 + 4 x 10 33 a; 21 , where V is denned to be 
dimensionless and of order unity, V = V/(\m pl ). This potential produces the fluctuation spectrum shown in Fig. 4 and the 
distribution of fluctuations shown in Fig. 5. 



FIG. 4. The spectrum of density fluctuations at horizon crossing, Sh, associated with the plateau model (70) as a function 
of the logarithm of the length scale L in units of pc. The plateau region seen in Fig. 3 produces the rapid rise in power, 
corresponding to ~ 10 32 g PBH production. 



FIG. 5. The solid line is the calculated final distribution of fluctuations (82) associated with PBH production from the 
plateau potential (70), Fig. 3. The dashed line is a stochastic numerical simulation of the the same model which consists of 
4 x 10 4 points. The calculated distribution is normalized to to the number of points and bin size of the numerical result. 



FIG. 6. A comparison of the calculated distribution of fluctuations for the plateau potential (70) with two Gaussian coun- 
terparts. The solid line is the analytic result (82). The dashed line is a Gaussian with the same mean and standard deviation 
as the calculated distribution. The short-dashed (dotted) line is a Gaussian with the same width at half maximum as the peak 
of the calculated distribution. Both Gaussians over-populate large fluctuations. 



FIG. 7. The same curves as in Fig. 6, but on a log scale to emphasize the error in the Gaussian assumption. Again, the solid 
line is the probability distribution (82) associated with the plateau potential (70). The two dashed lines are possible Gaussians 
that we may wish to compare our distribution to (which is which is not really important.) The distributions are plotted in 
terms of the number of standard deviations from the mean, where we have chosen the standard deviation of the peak of the 
distribution as the most conservative choice. A fluctuation of 6 — a would correspond to PBH production under the Gaussian 
assumption, but the actual distribution under-produces these fluctuations by ~ 10 150 , resulting in almost no PBH production. 



FIG. 8. The wiggle potential (83): V = 1+ [136.717]:r 3 — 0.05x. The bump in the path of the inflaton causes it to slow down, 
producing a spike in power shown in Fig. 9. 



FIG. 9. The spectrum of density fluctuations at horizon crossing, Sh, associated with the wiggle potential (83). The 
high-point in power corresponds to the production of ~ 10 28 g PBH's. 
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FIG. 10. The solid line shows the calculated distribution of fluctuations (94) associated with the wiggle model (83). The 
dashed line shows the results of our numerical simulation of this model, with 4 x 10 4 points. The distribution (94) was normalized 
to the simulation number and bin size, and is clearly non-Gaussian and skewed negative. 

FIG. 11. The cliff potential (95) which produces the blue spectrum shown in Fig. 12. The final distribution of fluctuations 
is nearly Gaussian since the slope of the potential has a sharp increase just after the spike in power (see Fig. 12). 

FIG. 12. The spectrum of density fluctuations at horizon crossing, 5h, associated with the cliff potential (95). The high-point 
in power corresponds to the production of ~ 10 20 g PBH's. 
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